DIPGEOPR_02

class: inverse,left, middle
background-image: url(data:image/png;base64,#background.png)
background-size: cover

##Módulo 2: Análisis espacial y Percepción Remota satelital óptica

### Análisis espacial

Matías Olea <br> <a href="https://orcid.org/0000-0003-0194-7784">
ORCID </a><br> 
matias.olea@pucv.cl</a><br>

.large[<b><a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Septiembre 2024</b>] <br>

---
class: center,middle
background-image: url(data:image/png;base64,#labgrs_logo.png)
background-size: 35%

---
## Contenidos

.pull-left[

1) Operaciones de consulta a un ráster.

2) Operaciones de transformación a un ráster.

3) Operaciones de alteración a un raster.

4) Operaciones de extracción de información.
  
  5) Análisis espacial entre ráster.
]

.pull-right[
<img src="data:image/png;base64,#https://raw.githubusercontent.com/allisonhorst/stats-illustrations/main/rstats-artwork/r_rollercoaster.png" width="650px"/>
]

---
## Modelo ráster

---
## Formatos ráster

Los modelos ráster pueden ser almacenados en múltiples formatos, los más comunes son:

--
- **GeoTIFF** (Geographical Tagged Image File Format). Fuente de datos en formato Tiff, con la particularidad de que es <u>capaz de almacenar información de referencia espacial</u>.

--
- **GMLJP2** (Geographic Markup Language JPEG2000). Muy parecido al formato JPEG pero con <u>mayores beneficios a la hora de comprensión de los datos</u> sin perder información o calidad.

--
- **NetCDF** (Network Common Data Form). Fuente de datos creada para el almacenamiento de información científica en varias dimensiones. Esta fuente de información se asocia principalmente a datos de tipos climáticos o productos satelitales. Su <u>gran capacidad de información</u> la hacen muy útil para compartir grandes volúmenes de información.

--
- **HDF** (Hierarchical Data Format). Al igual que el NetCDF, son fuentes creadas para almacenar <u>grandes volúmenes de información</u> científica espacial.

--
- **HGT** Formato asociado a modelo de elevación SRTM (Shuttle Radar Topography Mission Height) <u>caracterizado por almacenar valores enteros positivos y negativos de 2 bits</u>.
---
## Operaciones de consulta en R

``` r
library(terra)
elevacion_3672 <- rast("C:/TuDirectorio/")
plot(elevacion_3672)
```

---
## Operaciones de consulta en R

**Clase de objeto**

``` r
class(elevacion_3672)
```

```
## [1] "SpatRaster"
## attr(,"package")
## [1] "terra"
```

--
**Dimensiones**

``` r
c(nrow(elevacion_3672), ncol(elevacion_3672), nlyr(elevacion_3672))
```

```
## [1] 3601 3601    1
```

--
**Resolución**

``` r
res(elevacion_3672)
```

```
## [1] 0.0002777778 0.0002777778
```

---
## Operaciones de consulta en R

**Extensión**

``` r
ext(elevacion_3672)
```

```
##      xmin      xmax      ymin      ymax 
## -72.00014 -70.99986 -36.00014 -34.99986
```
---
## Operaciones de consulta en R

**Sistema de referencia y Coordenadas**

``` r
crs(elevacion_3672)
```

```
## [[1]]
##  [1] "GEOGCRS[\"WGS 84\","                                   
##  [2] "    DATUM[\"World Geodetic System 1984\","             
##  [3] "        ELLIPSOID[\"WGS 84\",6378137,298.257223563,"   
##  [4] "            LENGTHUNIT[\"metre\",1]]],"                
##  [5] "    PRIMEM[\"Greenwich\",0,"                           
##  [6] "        ANGLEUNIT[\"degree\",0.0174532925199433]],"    
##  [7] "    CS[ellipsoidal,2],"                                
##  [8] "        AXIS[\"geodetic latitude (Lat)\",north,"       
##  [9] "            ORDER[1],"                                 
## [10] "            ANGLEUNIT[\"degree\",0.0174532925199433]],"
## [11] "        AXIS[\"geodetic longitude (Lon)\",east,"       
## [12] "            ORDER[2],"                                 
## [13] "            ANGLEUNIT[\"degree\",0.0174532925199433]],"
## [14] "    ID[\"EPSG\",4326]]"
```

---
## Operaciones de consulta en R

**Tipo de dato**

``` r
datatype(elevacion_3672)
```

```
## [1] "INT2S"
```

**_Codificación_**

<div class="datatables html-widget html-fill-item" id="htmlwidget-ea77375aeefd43fb6da7" style="width:100%;height:auto;"></div>
<script type="application/json" data-for="htmlwidget-ea77375aeefd43fb6da7">{"x":{"filter":"none","vertical":false,"fillContainer":false,"data":[["1","2","3","4","5","6","7","8","9"],["LOG1S","INT1S","INT1U","INT2S","INT2U","INT4S","INT4U","FLT4S","FLT8S"],["FALSE (0)","-127","0","-32,767","0","-2,147,483,647","0","-3.4e+38","-1.7e+308"],["TRUE (1)","127","255","32,767","65,534","2,147,483,647","4,294,967,296","3.4e+38","1.7e+308"]],"container":"<table class=\"display\">\n  <thead>\n    <tr>\n      <th> <\/th>\n      <th>Datatype definition<\/th>\n      <th>Minimum possible value<\/th>\n      <th>Maximum possible value<\/th>\n    <\/tr>\n  <\/thead>\n<\/table>","options":{"pageLength":4,"autoWidth":true,"scrollY":true,"columnDefs":[{"orderable":false,"targets":0},{"name":" ","targets":0},{"name":"Datatype definition","targets":1},{"name":"Minimum possible value","targets":2},{"name":"Maximum possible value","targets":3}],"order":[],"orderClasses":false,"lengthMenu":[4,10,25,50,100]}},"evals":[],"jsHooks":[]}</script>
---

---

---

## Operaciones de consulta en R

Ahora, esta información que consultamos una por una, la podemos obtener de forma fácil y resumida solo imprimiento nuestro objeto:

``` r
elevacion_3672
```

```
## class       : SpatRaster 
## dimensions  : 3601, 3601, 1  (nrow, ncol, nlyr)
## resolution  : 0.0002777778, 0.0002777778  (x, y)
## extent      : -72.00014, -70.99986, -36.00014, -34.99986  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : s36w072.hgt 
## name        : s36w072
```

---

## Operaciones de consulta en R

**Consulta booleana**

``` r
consulta_250 <- elevacion_3672 
values(consulta_250) <- as.numeric(values(consulta_250) <= 250)
plot(consulta_250)
```

---
## Operaciones de alteración
** Uniones y Mosaicos**

Las uniones y mosaicos nos permiten unir más de un geodato ráster. Por principio, la union/merge nos permite unir cuando no hay zonas de traslape entre los ráster, y mosaico une cuando hay zonas de traslape. Teniendo en cosideración ello, la función **mosaic()** de terra, es más rápida que **merge()**.

---

## Operaciones de alteración
** Uniones y Mosaicos**

``` r
coleccion_elev <- list.files("TuDirectorio/", glob2rx("*hgt"), full.names=T) %>% 
  sprc()
mosaico_elev <- mosaic(coleccion_elev)
plot(mosaico_elev)
```

---
## Operaciones de alteración
**Recorte, Clip o Crop**

El recorte signica sustraer parte de la matriz ráster según una extensión definida. Podemos crear una extensión ingresando las <u>coordenas límite</u> (xmin, xmax, ymin, ymax), usando la <u>extensión de un geodato</u>, o <u>dibujando una extensión</u> sobre un plot usando la función **draw()**.

---
## Operaciones de alteración
**Recorte, Clip o Crop**

``` r
cuenca_maule <- vect("TuDirectorio/Cuenca_Rio_Maule.gpkg") %>% project(crs(mosaico_elev))
recorte_elev <- mosaico_elev %>% crop(ext(cuenca_maule))
plot(recorte_elev)
plot(cuenca_maule, add=T)
```

---
## Operaciones de alteración
**Máscara o Mask**

La aplicación de una máscara o enmascaramiento implica el sustraer los valores de un ráster dentro de los límites de un geodato, reclasificanto los pixeles fuera de estos límites como "No Dato", "not available" o "NA". Se aconseja recortar y enmascarar al área de interes.

---
## Operaciones de alteración
**Máscara o Mask**

Un elemento importante a considerar es hacerse la pregunta, ¿Qué ocurre cuando el límite de mi área de interés toca parte del pixel pero no lo cubre por completo? ¿lo incluye o no?

Existe un parámetro de tipo lógico dentro la función **mask()** llamado "touches". Cuando touches es **T** si toca el pixel, queda incluido en la máscara, si touches es **F** el pixel debe estar contenido en un 100% dentro del área de interes para ser incluido.

---
## Operaciones de alteración
**Máscara o Mask**

``` r
mascara_completa <- mosaico_elev %>% mask(cuenca_maule)
mascara_crop <- recorte_elev %>% mask(cuenca_maule)
```

---
## Operaciones de alteración
**Reclasificación**

La reclasificación significa reasignar nuevos valores a las celdas del ráster. Por lo general se utiliza cuando queremos pasar de valores continuos a valores discretos para agrupar por rangos en categorías específicas.

---
## Operaciones de alteración
**Reclasificación**

``` r
# Opción 1
reclass_1 <- mascara_crop # guardamos en otro objeto por si nos equivocamos

reclass_1[reclass_1 <= 500] <- 0
reclass_1[reclass_1 > 3000] <- -4
reclass_1[reclass_1 > 2000] <- -3
reclass_1[reclass_1 > 1000] <- -2
reclass_1[reclass_1 > 500] <- -1
```

---
## Operaciones de alteración
**Reclasificación**

``` r
# Opción 2
intervalos <- c(-Inf,500,0,500,1000,1,1000,2000,2,2000,3000,3,3000,Inf,4)
matriz_reclass <- matrix(intervalos, nrow=5, ncol=3, byrow=T)
```

```
##      [,1] [,2] [,3]
## [1,] -Inf  500    0
## [2,]  500 1000    1
## [3,] 1000 2000    2
## [4,] 2000 3000    3
## [5,] 3000  Inf    4
```
---
## Operaciones de alteración
**Reclasificación**

``` r
# Opción 2
reclass_2 <- classify(mascara_crop, matriz_reclass)
```

---
## Operaciones de alteración

**Remuestreo**

El remuestreo es una técnica que se utiliza para transformar la resolución de nuestro dato. Cuando queremos reducir el tamaño del pixel original, estamos hablando de **downscaling**, y cuando queremos aumentar el tamaño, **upscaling**.

Comunmente recurrimos al remuestreo cuando tenemos fuentes de <u>información con diferente resolución nativa</u>. Por ejemplo la información climática grillada del CR2 (CR2MET) tiene un pixel de 5 km. El modelo de elevación que estamos usando tiene una resolución de 30m. Las imagenes del sensor MODIS tienen una resolución de 250 m en algunos productos.

Por otra parte, tambien utilizamos el remuestreo cuando se realizan reproyecciones. Tópico que veremos en las operaciones de transformación.

Entre los métodos de remuestro más comunes están: Vecino más cercano (Nearest neighbor), Interpolación bilineal (Bilinear), Convolución cúbica (Cubic convolution) y por mayoría (Majority).

---
## Operaciones de alteración

**Remuestreo**

.pull-left[
**Vecino mas cercano (Nearest Neighbor)**

Este método asiga como valor de salida, el valor de la celda original cuyo centro esté posicionado en la nueva celda. Este método se utiliza cuando no queremos alterar el valor original de entrada sino que solo reposicionarlo, es el caso de las variables discretas, como por ejemplo los usos de suelo.

]

.pull-right[
<center><img src="data:image/png;base64,#vecino.png" width="500px"/></center>
]

---
## Operaciones de alteración

**Remuestreo**

.pull-left[
**Interpolacion bilinear (Bilinear)**

Este método no solo toma en cuenta el pixel de origen que ocupa el centro del nuevo pixel, sino que toma en consideración otros 3 y promedia los valores en porcentajes equivalentes a la superficie que cubren del nuevo pixel. Este método tiene a suavizar las transiciones de valores, se utiliza a menudo cuando hay ruido en los datos de entrada.

]

.pull-right[
<center><img src="data:image/png;base64,#BILINEAR.png" width="500px"/></center>
]

---
## Operaciones de alteración

**Remuestreo**

.pull-left[
**Convolución Cúbica (Cubic Convolution)**

Este método es parecido al bilineal pero en lugar de considerar 4 pixeles, considera 16. Es utilizado cuando existe exceso de ruido y tiene como resultado un suavizado mayor que el bilineal.

]

.pull-right[
<center><img src="data:image/png;base64,#CUBICA.png" width="500px"/></center>
]

---
## Operaciones de alteración

**Remuestreo**

.pull-left[
**Mayoria (Majority)**

Este método es muy similar al vecino más cercano, la gran diferencia es que en lugar de adquirir el valor del pixel que se posiciona en el centro del nuevo, adquiere el valor del pixel que cubre mayor parte de la superficie del nuevo.

]

.pull-right[
<center><img src="data:image/png;base64,#MAYORIA.png" width="500px"/></center>
]
---
## Operaciones de alteración

**Remuestreo**

<center><img src="data:image/png;base64,#plot_dem.png" width="1000px"/></center>
---
## Operaciones de alteración

**Remuestreo**

<center><img src="data:image/png;base64,#perfil.png" width="1000px"/></center>
---
## Operaciones de alteración

**Remuestreo**

---
## Operaciones de transformación

**Reproyección**

La reproyección implica reposicionar en el espacio XY nuestro geodato teniendo en consideración otro punto de referencia, proyección y/o coordenadas.

Si no conocemos el sistema de referencias y coordenadas, podemos revisar en [EPSG.IO](https://epsg.io/)

<center><img src="data:image/png;base64,#epsg.png" width="600px"/></center>
---
## Operaciones de transformación

**Reproyección**

.footnote[
<span style="font-size:9pt"> Fuente: T. Outtara et al (2004); Remote Sensing and Geosciences </span>
]

---
## Operaciones de transformación

**Reproyección**

``` r
reclass_utm <- project(reclass_2, "EPSG:32719", method="near", res=30)
crs(reclass_utm, describe=T)
```

```
##                    name authority  code
## 1 WGS 84 / UTM zone 19S      EPSG 32719
##                                                                                                                                           area
## 1 Between 72°W and 66°W, southern hemisphere between 80°S and equator, onshore and offshore. Argentina. Bolivia. Brazil. Chile. Colombia. Peru
##             extent
## 1 -72, -66, -80, 0
```

---
## Operaciones de transformación

**Cambio de tipo de dato**

``` r
datatype(elevacion_3672)
```

```
## [1] "INT2S"
```

``` r
summary(elevacion_3672)
```

```
##     s36w072      
##  Min.   :   3.0  
##  1st Qu.: 141.0  
##  Median : 240.0  
##  Mean   : 403.1  
##  3rd Qu.: 468.0  
##  Max.   :3031.0
```

``` r
writeRaster(elevacion_3672, filename="TuDirectorio/Elevacion.tif", datatype="INT1U")
```

---
## Operaciones de transformación

**Vectorización: poligonización**

``` r
reclass_pol <- as.polygons(reclass_utm, dissolve=T, values=T)
```

---
## Operaciones de transformación

**Vectorización: contornos o isolineas**

``` r
nivel <- c(250,500,1500,2500,3500)
elev_contorno <- as.contour(mascara_crop, levels=nivel) # creamos vector de lineas
plot(mascara_crop)
terra::contour(mascara_crop, levels=nivel, add=T) # nos permite plotear encima con etiqueta
```
<img src="data:image/png;base64,#DIPGEOPR_02_1_files/figure-html/unnamed-chunk-35-1.png" width="100%" />

---

## Operaciones de extracción de información.

**Cálculo de área**

``` r
frec_reclass <- freq(reclass_utm)
frec_reclass$area <- (frec_reclass$count * 900)/1000000
```

<div class="datatables html-widget html-fill-item" id="htmlwidget-9413fffadf4ea13d3cbd" style="width:100%;height:auto;"></div>
<script type="application/json" data-for="htmlwidget-9413fffadf4ea13d3cbd">{"x":{"filter":"none","vertical":false,"fillContainer":false,"data":[["1","2","3","4","5"],[1,1,1,1,1],[0,1,2,3,4],[12945730,2196707,4346529,3800804,121077],[11651.157,1977.0363,3911.8761,3420.7236,108.9693]],"container":"<table class=\"display\">\n  <thead>\n    <tr>\n      <th> <\/th>\n      <th>layer<\/th>\n      <th>value<\/th>\n      <th>count<\/th>\n      <th>area<\/th>\n    <\/tr>\n  <\/thead>\n<\/table>","options":{"columnDefs":[{"className":"dt-right","targets":[1,2,3,4]},{"orderable":false,"targets":0},{"name":" ","targets":0},{"name":"layer","targets":1},{"name":"value","targets":2},{"name":"count","targets":3},{"name":"area","targets":4}],"order":[],"autoWidth":false,"orderClasses":false}},"evals":[],"jsHooks":[]}</script>

---
## Operaciones de extracción de información.

**Cálculo de área**

.pull-left[
Abra el Landcover de la cuenca del Maule [(Zhao et al., 2016)](https://doi.org/10.1016/j.rse.2016.05.016), y calcule el área para cada una de las coberturas de uso de suelo descritas.

Utilice como referencia la tabla adjunta que detalla el significado de los valores del producto.

- ¿Cuál es la clase con mayor superficie?
- ¿Cuál es la superficie más pequeña en km2?
]

.pull-right[
<center><img src="data:image/png;base64,#landcover.png" width="350px"/></center>
]

---
## Operaciones de extracción de información.

**Cálculo de área**

``` r
landcover_maule <- rast("TuDirectorio/Landcover_Cuenca_Maule.tif")
frec_landcover <- freq(landcover_maule)
frec_landcover$value[which.max(frec_landcover$count)] # valor 410
min(frec_landcover$count)*900/1000000 # 0.0225 km2
```

```
## [1] 410
```

```
## [1] 0.0225
```

---
## Operaciones de extracción de información.

**Extracción de valores de los pixeles**

``` r
hist(mascara_crop, main="Histograma elevación cuenca Maule")
```

---
## Operaciones de extracción de información.

**Extracción de valores de los pixeles**

``` r
valores_elev <- values(mascara_crop) %>% na.omit # omitimos las celdas NA
mean(valores_elev)
```

```
## [1] 875.3846
```

``` r
median(valores_elev)
```

```
## [1] 388
```

---
## Operaciones de extracción de información.

**Estadística zonal**

Para este ejercicio extraeremos el valor promedio de elevación por ciudad de la provincia de Talca de la región del Maule.

``` r
comunas_maule <- vect("TuDirectorio/Region_Maule_Difrol_latlong.gpkg")
provincia_talca <- subset(ciudades_maule, ciudades_maule$NOM_PROVIN == "TALCA")
media_elev_comunas <- terra::extract(mosaico_elev, provincia_talca, 
                                     fun=mean, method="simple", touches=T)
media_elev_comunas$ciudad <- provincia_talca$URBANO
```

---
## Operaciones de extracción de información.

**Estadística zonal**

<div class="datatables html-widget html-fill-item" id="htmlwidget-dd94b61acaf16fb16919" style="width:100%;height:auto;"></div>
<script type="application/json" data-for="htmlwidget-dd94b61acaf16fb16919">{"x":{"filter":"none","vertical":false,"fillContainer":false,"data":[["1","2","3","4","5"],[1,2,3,4,5],[45.98044032882925,200.7756241441062,115.2182389937107,102.6044285692196,96.59722397476341],["CONSTITUCIÓN","SAN CLEMENTE","MAULE","TALCA","CULENAR"]],"container":"<table class=\"display\">\n  <thead>\n    <tr>\n      <th> <\/th>\n      <th>ID<\/th>\n      <th>s35w071<\/th>\n      <th>ciudad<\/th>\n    <\/tr>\n  <\/thead>\n<\/table>","options":{"columnDefs":[{"className":"dt-right","targets":[1,2]},{"orderable":false,"targets":0},{"name":" ","targets":0},{"name":"ID","targets":1},{"name":"s35w071","targets":2},{"name":"ciudad","targets":3}],"order":[],"autoWidth":false,"orderClasses":false}},"evals":[],"jsHooks":[]}</script>

---
## Análisis espacial entre ráster.

**Algebra de mapas o Algebra ráster**

---

## Análisis espacial entre ráster.

**Algebra de mapas o Algebra ráster**

En un caso hipotetico donde haya que definir sitios de interés cientifico o para la conservación, se definieron las variables de EVI (Indice de Vegetación mejorada), la elevación y el Landcover de la cuenca del río Maule para generar zonas de interés.

Teniendo en consideración un puntaje de 0 para zonas menos aptas, y 4 para más aptas, reclasifique sus variables de acuerdo a los siguientes criterios:

Una vez reclasificado, multiplique las 3 variables y revise cuales son las áreas con mayor puntaje.

---

## Análisis espacial entre ráster.

**Algebra de mapas o Algebra ráster**

``` r
## Apertura de Datos
elev_maule <- rast("TuDirectorio/ELEV_CUENCA_MAULE.tif")
land_maule <- rast("TuDirectorio/Landcover_Cuenca_Maule.tif")
evi_maule <- rast("TuDirectorio/EVI_CUENCA_MAULE.tif")

## Matrices de reclasificación
class_elev <- c(-Inf,0,0,0,1000,1,1000,2000,2,2000,3000,3,3000,Inf,4)
matriz_elev <- matrix(class_elev, nrow=5, ncol=3, byrow=T)

class_evi <- c(-Inf,0,0,0,0.15,1,0.15,0.3,2,0.3,0.45,3,0.45,Inf,4)
matriz_evi <- matrix(class_evi, nrow=5, ncol=3, byrow=T)

class_land <- c(0,150,0,200,239,4,240,260,1,300,390,2,400,480,4,500,700,3,799,Inf,0)
matriz_land <- matrix(class_land, nrow=7, ncol=3, byrow=T)
```

---

## Análisis espacial entre ráster.

**Algebra de mapas o Algebra ráster**

---

## Análisis espacial entre ráster.

**Algebra de mapas o Algebra ráster**

.pull-left[

Boxplot de los pixeles

``` r
boxplot(Area_interes,horizontal=F)
```

]

.pull-right[
Cálculo del percentil 95

``` r
quantile(values(Area_interes) %>% na.omit(),
         probs = 0.95)
```

```
## 95% 
##  24
```
]

---
## Análisis espacial entre ráster.

**Algebra de mapas o Algebra ráster**

---

### Bibliografía

2022 Jerry Davis. Introduction to Environmental Data Science. Capitulo 8.Raster Spatial Analysis. https://bookdown.org/igisc/EnvDataSci/raster-spatial-analysis.html

2022 Robert J. Hijmans. terra.pdf. https://cran.r-project.org/web/packages/terra/terra.pdf

2012 Essentials of Geographic Information Systems. Chapter 8.Geospatial Analysis II: Raster Data. https://saylordotorg.github.io/text_essentials-of-geographic-information-systems/s12-geospatial-analysis-ii-raster-.html

---
class: inverse middle
![logo](data:image/png;base64,#LABGRS-PUCV.png)